//
//  main.cpp
//  EMT_V2_Crisis
//
//  Created by Duong Ngo on 6/5/17.
//  Copyright © 2017 Duong Ngo. All rights reserved.
//

#include <iostream>
#include "global.hpp"
#include "mapping.hpp"
#include "read_steady_state.hpp"
#include "foc_problem.hpp"

int main()
{
    // Global mapping and vector showing position
    find_original_map(original_map);
    find_equation_map(eqn_map);
    find_positions_of_state_and_forward_variables(state_position, forward_position, var_role);

    // Read steady state
    std::vector<double> ss(var_foc);     //Steady State
    read_steady(ss);
    print_vector(ss);

    
    
    // Initial guess
    std::vector<double> x_init((T+1)*var_foc);
    initial_guess_as_steady_state(ss, x_init);
    
    // Starting state
    starting_state[0] = ss[original_map.at("s.n")];
    starting_state[1] = ss[original_map.at("s.m")];
    starting_state[2] = ss[original_map.at("s.bh")];
    starting_state[3] = ss[original_map.at("s.k")];
    starting_state[4] = ss[original_map.at("s.Rm")];
    
    // Exogenous shock
    std::vector<double> kappa_start(T+1), kappa_end(T+1);
    for (int t=0; t<T+1; ++t){
        kappa_start[t]= kappa;
    }
    vec_kappa = kappa_start;
    
    kappa_end[0] = kappa_shock;
    for (int t=1; t<T+1; ++t){
        kappa_end[t]= rho*kappa_end[t-1]+ (1-rho)*kappa;
    }
    
    // Solve a sample
    foc_solve(x_init, ss);
    

    // Homotopy
    int N_iter = 100;
    int dem =0;
    double chi=0; //Homotopy coefficient
    double acceptable_error = 1e-6;
 
    while (dem <= N_iter && max_error < acceptable_error){
        printf("*********************The iteration: %2u *********************\n", dem);
        chi= double(dem)/N_iter;
        for (int t=0; t<T+1; ++t){
            vec_kappa[t]= chi*kappa_end[t]+ (1-chi)*kappa;
        }
        foc_solve(x_init, ss);
        if (max_error < acceptable_error){
            x_init = result;
            store_output(result, chi);
        }
        else {
            for (int j=1; j<= N_iter; ++j){
                chi = double(dem-1)/N_iter+ double(j)/N_iter/N_iter;
                for (int t=0; t<T+1; ++t){
                    vec_kappa[t]= chi*kappa_end[t]+ (1-chi)*kappa;
                }
                foc_solve(x_init, ss);
                if (max_error < acceptable_error){
                    x_init = result;
                    store_output(result, chi);
                }
                else {
                    break;
                }
            }
        }
        ++dem;
    }

 
    return 0;
}
